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Abstract. We study, by numerical simulations and semi-rigorous arguments, a two- 
parameter family of convex, two-dimensional billiard tables, generalizing the one-parameter 
class of oval billiards of Benettin-Strelcyn |BS781. We observe interesting dynamical 
phenomena when the billiard tables are continuously deformed from the integrable circular 
billiard to different versions of completely-chaotic stadia. In particular, we conjecture that 
a new class of ergodic billiard tables is obtained in certain regions of the two-dimensional 
parameter space, when the billiards are close to skewed stadia. We provide heuristic arguments 
supporting this conjecture, and give numerical confirmation using the powerful method of 
Lyapunov-weighted dynamics. 
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1. Introduction 

Billiard models are a class of Hamiltonian dynamical systems which exhibit the full range of 
behaviour from completely integrable to completely chaotic dynamics I CM061 . They consist 
of a point particle which collides elastically with the walls of a bounded region, the billiard 
table; the shape of the table determines the type of dynamics which is observed. 

Several classes of billiard tables have been studied which interpolate between completely 
integrable and completely chaotic dynamics, including both one-parameter |BS78 HW83] 
IBunOll |MKd96l |Fol02 1 and two-parameter IIDRW96 I families. These allow us to observe 
the transitions by which the typical phase space, which is a mixture of ergodic, chaotic 
components and regular KAM islands, evolves from one extreme behaviour to the other. 

In the context of two-dimensional billiards, a popular starting point of such investigations 
is the Bunimovich stadium, which consists of two semicircular arcs connected by two parallel 
segments. The case when these two segments are non-parallel (and, correspondingly, one 
of the arcs is shorter, with the other being longer than a semicircle) is often referred to as 
the skewed stadium or squash billiard table. Stadia are known to be ergodic and hyperbolic; 
however, as a consequence of so-called "quasi-integrable" phenomena, the hyperbolicity is 
non-uniform, and the dynamical behaviour is very sensitive to perturbations of the boundary. 
There is an abundance of literature on stadium billiards; some works that are relevant to our 
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discussion are refs. IIBun74l lMar98l ICZ05I IBG06I IDG09I : see also section [TT] for a more 
detailed description of stadia. 

In this paper, we study a two-parameter set of two-dimensional billiard tables which 
generalizes the one-parameter family of oval billiards studied in |BS78 HW83 1 in a particular 



way; see section 2.2 for an explicit description. Our models constitute a subcase of a rather 



general class of billiards introduced, but not studied in detail, by Hayli and Dumont in BHD86I . 
Our class is of particular interest since it includes as limiting cases an entire family of systems 
which are known to be ergodic, the skewed stadium billiards. Here these are generalised by 
deforming the sides of the stadia to circular arcs. 

As in previous works on billiards formed from piecewise smooth components IIBS78I 
IHW83I IHD86I IDRW96I ILMROll ILMR02I ILMPR06I iLSTOl . for a large set of parameter 
values, we find coexistence of stability islands and chaotic components in phase space. 
However, we also find numerically that billiards which are sufficiently close to the limiting 
skewed stadia appear to have no remaining stability islands - the phase space is filled by a 
single large chaotic, ergodic component. This motivates our Conjecture |3. 1 [ concerning a new 
class of ergodic billiard tables. Similar conjectures have been suggested for other billiard 
models of similar type - see, e.g., IIDRW96 1 . 

Our case is, however, different from previous studies in several respects. Firstly, the 
new class introduced in this paper and conjectured (for an open set of parameter values) to 
be ergodic consists of convex planar billiards. The issue of ergodicity versus KAM islands 
in convex billiards has been in the focus of continued interest for several decades. On the 
one hand, by Lazutkin's observation in flLaz731 . and its strengthening according to Douady 
nDou82l . a convex planar billiard with at least boundary cannot be ergodic. Furthermore, 
recent results by Bunimovich and Grigo IIGri091lBG1 0l show that elliptic islands arise in 
stadium-like billiards (billiard tables constructed from stadia by replacing the discontinuity 
of the curvature with a smoothening). On the other hand, several classes of convex 
planar billiards (with some discontinuity points of the curvature) are proved to be ergodic 
(see section [TT] for a partial list of references). A common feature of these examples is the 
defocusing mechanism, which requires that whenever a narrow beam of (initially) parallel 
rays completes a series of consecutive reflections on one of the smooth focusing components 
of the billiard boundary, it must pass through a conjugate point and become divergent before 
the next collision with the curved (non-flat) part of the boundary; see [ICM06I IGri09l for a 
more detailed description. 

It is easy to see that in our examples defocusing cannot take place in this sense. Hence 
the mechanism that produces the (conjectured) ergodic behaviour must be different. We are 
aware of two other examples of ergodic planar billiards with focusing boundary components 
where defocusing is violated |BL08 BdM09|; however, non-convexity of the billiard domain 
plays an important role in both cases. Similarly, numerical studies in IIDRW96IIDB01II suggest 
ergodicity only for certain «o«-convex domains. 

Secondly, even though a rigorous proof is currently not available, we give, in addition 
to the simulated phase portraits, further evidence which strongly supports our conjecture. 



Heuristic arguments are provided in section 4.1 which rely on the similarity of the 
dynamics with those of the skewed stadia, and on the explicit analysis of sliding trajectories. 
Furthermore, the absence of islands is also tested numerically by the powerful method of 



Lyapunov-weighted dynamics in section 4.2 



The rest of this paper is organized as follows. In section|2] we define the two-parameter 
family of generalised squashes and highlight its differences from related models. Section [3] 
reports numerical results concerning the dependence of the dynamics on the two geometrical 
parameters, and presents the main conjecture on the existence of a new class of ergodic 
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billiards. Section |4]presents heuristic arguments and further numerical evidence supporting 
the conjecture; and in section[5]we draw some conclusions. 

2. The model 

We begin by defining the class of generalised squash billiard models. 
2.1. Convex billiard tables 

Consider a convex, compact domain Q <zM} bounded by a closed, piecewise-smooth curve 
T = dQ. The motion of a point particle that travels along straight lines with constant speed 
in the interior of Q, and bounces off elastically (angle of reflection equals angle of incidence) 
when reaching the boundary F, is referred to as billiard dynamics. 

We investigate these dynamics in discrete time, that is, from collision to collision. The 
phase space is then the cylinder M = F x [0, tt], with M 3 x = {k,<p), where the configurational 
coordinate is the arc length k along the boundary, which satisfies < < |F| and describes the 
point of the closed curve F at which the collision takes place, while the velocity coordinate 
< <p < TT describes the angle that the outgoing velocity makes with the (positively oriented) 
tangent line to F at the point k. 

Given x e M, the position and velocity at the next collision are uniquely determined, so 
that the billiard map T : M M is well-defined. It is usual to visualize M as a rectangular 
domain in the plane, and the consecutive points along a trajectory of T as points in this 
domain. T has a natural invariant measure pL, which is absolutely continuous with respect 
to Lebesgue measure on M, given by 

dpi = const, sin ^dkd^. 

For further material on billiards in general we refer to the monographs IICM06llTab951 . 

The billiard map may show a surprisingly wide variety of dynamical phenomena for 
different choices of the billiard table F. The best known case is the billiard in a circle: for this 
geometry the dynamics are integrable: the angle of incidence (p is an integral of motion, the 
values of which label the invariant curves l,CM06l . 

If F is a C^^^-smooth closed curve, then for trajectories in the vicinity of the boundary the 
dynamics resemble, to some extent, those of the circular billiard: Lazutkin fLaz73l showed 
that a positive measure set of the phase space in a neighborhood of the boundaries ^ = and 
(p = ;r is foliated by invariant curves. Later, Douady ( l|Dou82| ) lowered the requirement for 
Lazutkin's result to C^-smooth boundaries. 

However, if F has less smoothness, then the billiard can be completely chaotic and 
ergodic. The first examples of such billiard tables are the celebrated stadia: 

• the straight stadium is formed by two identical semicircles, joined at their endpoints by 
two parallel lines along their common tangents; 

• the skewed stadium, or squash table, is formed by two circular arcs of different radii 
r < R, joined at their endpoints by non-parallel straight lines along their common 
tangents. 

Stadia were introduced and their chaotic character first studied by Bunimovich in his 
famous paper |Bun74|. It is known that for these tables the billiard map T is completely 
hyperbolic, i.e. there is one strictly positive and one strictly negative Lyapunov exponent for 
;U-a.e. X e M, and T is ergodic with respect to jU; see ICZ05I ICM06I ICZ071 for a detailed 
description of stadia. 
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Stadia have three important characteristic features: 

• The table boundary Y is only piecewise smooth: at the intersection points of the 
circular arcs and straight lines, the curvature of the boundary (the second derivative) 
is discontinuous. The resulting singularities play a crucial role in the dynamics of the 
billiard map, both for the stadia and for the tables investigated below. 

• The map T has many periodic points, all of which are hyperbolic. In the case of 
the circular billiard, all periodic points are parabolic. We will see below that in 
the generalised squash tables, the third class, elliptic periodic points, may also arise, 
typically giving rise to KAM islands in their vicinity. 

• The billiard map T in stadia is only non-uniformly hyperbolic. The reason for this is the 
presence of so-called quasi-integrable - sliding, bouncing and diametrical - trajectory 
segments of arbitrary length. When the geometry is perturbed, these quasi-integrable 
phenomena may (or may not) create islands of integrability; see below for further 
discussion. 

The mechanism that is responsible for the chaotic behaviour in stadium-shaped tables, 
known as the defocusing mechanism, has also been observed and studied in other classes 
of billiard tables in works of Bunimovich, Donnay, Markarian, Szasz and Wojtkowski 
| |Woj86| IMar88l IDon91l IBun92l ISza92l IMar93L However, the geometry of the billiards 
studied in this paper is quite different from the geometry of any of these classes. 

2.2. The two-parameter family of generalised squashes 

The family of convex billiard tables studied in this paper is described by two parameters, b 
and c. The table is built on a trapezoid; on each side of the trapezoid is placed a circular arc 
joining its end-points, and adjacent arcs are constrained to meet with common tangents. The 
parameter < b < 1 specifies the geometry of the trapezoid, while 1 < c < oo specifies the 
ratio of the radii of the arcs. 

More precisely, the base and the two sides of the trapezoid are fixed to have unit length, 
and b is the length of the top, which is parallel to the base. The two extreme cases b — and 
b — 1 correspond to the equilateral triangle and the square, respectively; see figure[TJa). 

The table is left-right symmetric, so that the circular arcs corresponding to the two sides 
have the same radius. The radii of the arcs of the base, the sides and the top are denoted R, 
Roo and r, respectively. As adjacent arcs are required to have a common tangent line where 
they meet, the shape of the table is determined by the value of any one of the radii - i.e., 
for a given trapezoid (value of b), the geometry can be parametrized by a single quantity. 
For convenience, we choose this parameter as the ratio c :— Ro<,/R; see figure [T| For brevity, 
throughout the paper the arcs corresponding to the base, the sides and the top are referred to 
as the bottom arc, the "almost-flat" arcs, and the top arc, respectively. The construction of 
the tables is detailed in the Appendix; they can be viewed as a subfamily of a general class of 
billiards introduced in IIHD86I . 

Note that the case c = 1 corresponds to a circle for any value of b, the case b = I and 
c = °° corresponds to the straight stadium, while b < I and c — °a gives a skewed stadium, 
with the amount of skewness depending on the value of b. Thus, by changing c, we may 
continuously deform the integrable dynamics of the circular billiard into a completely-chaotic 
billiard. 

In this paper we concentrate on the role of the parameter b in this transition, which 
determines the relative width of the top and bottom sections of the billiard. The transition for 
the case b — 1 (oval billiards based on a square) has already been studied in the literature. 
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(a) Trapezoids for b = 0, 0.2, (b) Effect of the parameter c, for b = 0.4. 
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Figure 1. Variation of tlie geometry of tiie generalised squash billiards as a function of the two 
parameters b and c. 

although with a different parametrisation IIBS78I IHW83I . where conclusions were drawn 
based on the phase portraits obtained. The picture for b = \ can be summarized as follows. 
As soon as c is finite - that is, the table differs from the straight stadium - ergodicity of 
the billiard map is ruined. More precisely, elliptic islands - regions of positive Lebesgue 
measure, foliated by invariant curves and concentrated around elliptic periodic points - appear 
in the phase space. For intermediate values 1 < c < oo, coexistence of such elliptic islands 
and ergodic components of positive measure is observed. For large enough c, there is just a 
single such chaotic component, which still coexists with elliptic islands, while the number of 
chaotic components increases once c is decreased below certain threshholds. The references 
IIBS78||HW831 concentrated on this phenomenon of splitting of chaotic components. 

In the present paper our main goal is different: by extending our investigation lob < 1, 
we observe distinct phenomena for finite, but large enough, c. 

2.3. Quasi-integrable phenomena 

Certain differences between the b = \ and b < \ cases involve quasi-integrable trajectories, 
which, as mentioned above, make stadia only non-uniformly hyperbolic. These quasi- 
integrable phenomena have different geometric origins for straight and skewed stadia. Sliding 
- where the trajectory collides almost tangentially at a circular arc for many consecutive 
occasions - is present in both of these tables; however, bouncing - where the trajectory 
travels back and forth between two flat boundary components - is characteristic only of the 
straight stadium; cf. figure |2];a). More precisely, the time of bouncing is uniformly bounded 
for skewed stadia, whereas it can last arbitrarily long in the straight case. When c <oo^ the 
geometry is perturbed qualitatively at the almost-flat sides, which is exactly where bouncing 
takes place. This makes the quasi-integrable bouncing phenomena integrable - an elliptic 
island is created - if b — 1. However, when b <\, the role of bouncing is negligible, and so 
the changes are not so pronounced. 

Finally, in skewed stadia a third type of quasi-integrable motion may also arise, in 
addition to bouncing and sliding. When a circular arc is longer than a semi -circle, trajectories 
may travel back and forth between two opposite points of the circular arc, with velocity almost 
parallel to the diameter, for an arbitrarily long time; cf. figure |2|b). However, for finite c, 
the qualitative changes in the geometry do not affect this diametrical motion, which thus 
remains a source of quasi-integrable phenomena, and does not create islands of stability. See 
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(a) Bouncing and sliding orbits. (b) Diametiical orbits. 

Figure 2. Different types of quasi-integrable pfienomena present in stadia. 

|[Mar98' ' CZ05I IBG06I ICZ08I for discussions of the effect of quasi-integrable phenomena on 
the ergodic properties of stadia. 

3. Parameter-dependence of dynamics 

In this section, we survey the types of dynamical behaviour observed in different regions of 
parameter space. The numerical methods that were used are briefly described in the Appendix. 

3.1. Parameter space 

We now proceed to investigate the dependence of the dynamics on the parameters of the 
system, by systematically scanning the two-dimensional parameter space and identifying the 
different phenomena observed. Numerical results for other billiard models of similar types, 
formed by piecewise-smooth curves, can be found in I BS78..HW83 , HD86 , DRW96 , L MROfl 
ILMR02I ILMPR061 iLSTOl . 

For not too small c (for c > 1 .2 in all cases), ergodic components of positive Lebesgue 
measure seem to fill most of the phase space. An initial condition chosen randomly in such a 
component has an orbit which fills the component densely. For c > 1.5, we observe a single, 
dominant ergodic component, whereas for smaller values of c several ergodic components 
coexist. For b—\ this "splitting of chaotic components", studied in |BS78, HW83|, occurs at 
c ~ 1.35; for < 1 we observe similar phenomena, although restricted to a smaller range of 
c. 

If the initial point is chosen in the complement of the large ergodic component(s), then we 
observe the typical mixed phase space characteristic of KAM theory: regions centered around 
elliptic periodic points where a positive measure set foliated by invariant curves and chaotic 
orbits coexist in a topologically complicated way. We refer to these regions as "islands". 

3.2. Stability regions of periodic orbits 

To gain understanding of the observed phenomena, it is useful to start from the known 
behaviour of the case b= I, i.e. the oval billiards studied in IIBS78IIHW83J . 
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Figure 3. The period-2 orbit and corresponding phase-space structures for c=100 with b=\ 
(left) and b = 0.99 (right). 



For these billiards, an important feature of the phase portrait is a period-2 orbit which is 
present for any finite value of c, and which bounces perpendicularly between the midpoints of 
the almost-flat sides. Its stability properties can easily be calculated (cf. Appendix A): when 
c < oo (that is, Roo < °°), it is elliptic, and gives rise to a KAM island. The island around this 
orbit remains fully extended in the k coordinate along the almost-flat arcs, and it shrinks only 
in the (p direction as we increase the value of c; see figure|3ja). 

For b < 1, however, the island shrinks both in the k and in the (p directions as c is 
increased (see figure [3|b)), and finally disappears for large enough c. The reason for this is 
that the almost-flat arcs are now placed along two non-parallel sides of the trapezoid, so that 
for large enough c there no longer exist any points on these two arcs which have mutually 
parallel tangent lines, and hence there is no period-2 orbit. In this case, the bottom arc is 
forced to be longer than a semi-circle, due to the orientation of the tangent lines at the points 
where the bottom and almost-flat arcs join. The limit of existence of this period-2 orbit is thus 
when the bottom arc is exactly a semi -circle, which occurs at c = co{b) :— 

For c > co{b), the period-2 orbit and its corresponding island disappear, but other islands, 
corresponding to periodic orbits of higher periods, appear in certain regions of the {b,c) 
plane. These islands generally appear in a certain window of values of c for a given, fixed 
b. Such stability regions for certain orbits of low period are shown in figure |4] For example, 
in the (blue) B region in the parameter space of figure |4] there exists an island around the 
period-4 orbit shown at the top right, for which the bounces on the bottom and left arcs are 
perpendicular (see also figure[5ja)). 



3.3. Geometric destruction and period doubling of periodic orbits 

We proceed to investigate some features of the phase portrait, in particular, the regions of 
existence and stability for the periodic orbits. Bifurcation phenomena in piecewise-smooth 
systems are the subject of intensive studies - see e.g. IIDBBCK08I ISM09I and references 
therein. It is beyond the scope of this paper, and in any case not our main goal here, to 
provide a detailed study of the bifurcations occurring in the system; instead, we give a short 
description of the main features we have observed (see BDRW96 1 for a similar discussion). 

Based on our observations, there are two kinds of dynamical phenomena that mainly 
determine whether an island of stability is formed around a periodic orbit. On the one hand, 
the singularities of the boundary play an important role, since they influence the very existence 
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parameter b 

Figure 4. Regions of stability of certain types of periodic orbit. Each shaded region depicts the 
numerically-determined region of parameter space in which the corresponding labelled orbit 
type (with the same topology, or coding) is stable. 



of periodic orbits. To describe these orbits succinctly, we introduce a coding of four letters 
{Z7,r,f,/}, denoting bounces on the bottom, right, top and left arcs, respectively, and describe 
periodic orbits according to a finite code of length equal to the period of the orbit. Even 
though this code is not unique, it is useful since it highlights the role of the singularities, as 
far as the existence and stability features of periodic orbits are concerned. 

For example, consider again the period-2 orbit with code (r, /), which consists of 
consecutive perpendicular bounces on the two almost- flat arcs, and which exists for c < co{b). 
As c tends to co{b) from below, the locations of the bounces on the sides move closer 
and closer to the ends of the almost-flat arcs, and finally reach the lower comer points 
(singularities) when c = co{b). For larger values of c, this period-2 orbit with code {r,l) 
thus ceases to exist, since there are no points on the almost flat arcs with mutually parallel 
tangent lines. We refer to this as geometric destruction of the periodic orbit. 

However, geometric destruction is not the only mechanism by which an island 
corresponding to a periodic orbit may disappear For example, figure |5]^a) shows a period- 
4 orbit, around which an island of stability exists in the range 6.7 < c < 8.7 for for b — 0.4. 
When c ~ 8.7, the orbit loses stability, even though its points of collision are still located 
far from the singularities. In this case, the disappearance of the island is rather related to 
the stability properties of the periodic orbit. KAM theory can only account for islands near 



elliptic periodic points, i.e. if the stability parameter (see equation (B.l i and the following 
paragraph) satisfies \s\ < 2. 

In the smooth setting, it is known that when a periodic orbit loses stability (\s\ increases 
above 2) period-doubling bifurcations may be observed; see eg. IIMac93ll . As observed in 
|[bRW96|, this can also occur in billiard-type dynamics, provided the periodic orbit stays 
away from the singularities. 

Consider again the period-4 orbit studied above for b — 0.4. As c increases from c ~ 6.7 
to c ~ 8.7, we find that s decreases from 2 to —2; see figure |7] When c reaches the critical 
value 8.7, the shape of the island changes and it splits into two components, the centers of 
which are consecutive points of a period-8 orbit, shown in figure |5|b). The changes in shape 
of the corresponding islands which surround the periodic orbits is shown in figure [6] 



Chaos and stability in a two-parameter family of convex billiard tables 



9 




(a) Period-4 orbit with c = 8.59. (b) Peiiod-8 orbit with c = 8.61. 

Figure 5. Period-doubling bifurcation of a period-4 orbit for b = 0.4 as c is varied. The (red) 
thin straight lines indicate the singularities in the boundary of the squash, where the circular 
arcs are joined. The arrows indicate the direction in which the trajectory moves as c increases. 
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Figure 6. Phase space structures for the peiiod-doubling bifurcation of the elliptic periodic 
orbit shown in figure |5j a), highlighting the splitting of the corresponding island. Parameter 
values b = 0.4 and c = 8.59, 8.61, 8.63 and 8.69 from left to right. 

The period-8 orbit is elliptic even for larger c, when the original period-4 orbit is already 
hyperbolic, so the island is formed around this elliptic orbit of double period for these 
parameter values. Similar period-doubling bifurcations occur for other periodic orbits whose 
disappearance is not due to geometric destruction. 

As c is increased even further, the geometry of the period-8 orbit is deformed, and 
soon the location, for at least one of the bounces, reaches a corner, resulting in geometric 
destruction of the island. 

Summarizing, we can see that islands are subject to the combination of two effects: 
geometric destruction caused by the singularities, combined with the generic features 
characteristic of smooth systems with mixed phase space. Indeed, the presence of the 
singularities modifies the picture significantly. See also IIDRW96I IAMS09I on similar 
phenomena in the billiard setting. 
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Figure 7. Dependence of the stability on c for the period-4 and period-8 orbits shown in 
figure|5] 

3.4. Main conjecture: a class of ergodic convex billiards 

In this section we concentrate on the case of large c for some b < \ fixed, and formulate our 
main conjecture concerning this region of parameter space. 

The tendency is that as c increases, the islands of stability are formulated around orbits of 
higher period. Typically, orbits of higher period are more sensitive to geometric destruction, 
since the presence of more bounces leads to more possibilities for one of the collision points 
to collide with a singularity when the parameters are changed. Thus the area in phase space of 
higher-period islands tends to decrease with increasing period. For this reason, the region in 
the parameter space where period-doubling bifurcation phenomena can actually be observed 
is quite limited. 

Furthermore, as the parameter c is increased even further, it appears that all islands 
disappear, and that the system thus becomes ergodic. These numerical observations motivate 
the following conjecture: 

Conjecture 3.1 For any b < \, there exists c{b) such that for any pair of parameters {b,c) 
with c > c{b) the billiard dynamics are ergodic. 

That is, once the generalised squash billiards are close enough to skewed stadia, we conjecture 
that they are ergodic. There is an obvious lower bound c{b) > cq [b) on the region of ergodicity 
in parameter space; both quantities tend to oo as tends to 1. 

Similar conjectures have been made for certain regions of parameter space in other 
billiards with piecewise-smooth boundaries |DRW96|. Here, however, we provide heuristic 
arguments, and confirmation using a more powerful numerical method, in the next section. 

4. Heuristic and numerical evidence for Conjecture |3.1| 

In this section, we support Conjecture |3 . 1 1 with heuristic and further numerical arguments. 
Even though a rigorous proof, which appears to be a technically challenging task, is currently 
not available, we believe that we capture the key phenomena. We also apply the powerful 
numerical method of Lyapunov-weighted dynamics to perform a more stringent search for 
elliptic islands. 
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4.1. Heuristic arguments 

We compare the dynamics of the model for large c with the limiting case c = oo, i.e. the skewed 
stadia, which has been studied extensively in the literature, for any value of b, and is known 
to be ergodic; see, for example, BCZ05L 



Comparison with the limiting case c ^ °o. Consider the subset M C M consisting of phase 
points X = {k, (p) for which (i) k belongs to one of the circular arcs (not to the almost-flat 
ones), and (ii) the (p coordinate ensures that for {k' ,(p') = T{k,(p), the point k' belongs to an 
arc different from the one on which k is located. It is proved in 0CZO5I that for c = °o, the first 
return map to M is uniformly hyperbolic. 

Now, for fixed b the billiard tables corresponding to finite c = c ^ 1 and to c = oo are 
(piecewise) C^-close. Let us denote the corresponding billiard maps by Tg and Toe, and the 
return maps to M by and T„, respectively. Now we recall some properties of the billiard 



map from I.CM061 ; see also formula ( |B.1| in the Appendix. C^-closeness of the billiard 
tables would imply that the maps Tc and Too are -close; however, the billiard map may have 
unbounded derivatives corresponding to tangential collisions. Thus, and Too are C^-close 
unless ^ ~ or ^ ~ TT. 

We would like to conclude that (apart from tangencies) Tc is also -close to Too, and thus, 
Tc is also uniformly hyperbolic (as uniform hyperbolicity is a C'-open property). However, 
closeness of the T's implies closeness of the 7"s only if phase points x EM are considered for 
which Tx = T"^^'^x such that n{x), the number of iterates needed to return to M, is uniformly 
bounded. Thus, we need to consider the complement of this set: points with unbounded return 
time, which are exactly the quasi-integrable trajectories. At this point the differences between 
the b < I and b ~ 1 cases play an important role, as follows. 

(i) For bouncing points x — {k, (p) G M, returns to M consist of a high number of consecutive 
collisions on the almost flat arcs. If b = 1, then such orbits may spend an unbounded 
number of iterations bouncing close-to-perpendicularly on the almost-flat arcs. However, 
if b < 1, then the number of such bounces, and thus the time needed to return to M, is 
uniformly bounded (the bound, of course, depends on the actual value of b, which affects 



the value of c{b) in Conjecture 3.1 



(ii) Diametrical quasi-integrable motion, when the trajectory moves back and forth between 
two opposite points of a circle (see the end of section |2]l, may also correspond to 
unbounded return time. This phenomenon, however, is restricted to the bottom arc of the 
b < I case, which has the same geometry for Tc and Too. More precisely, the almost-flat 
arcs play a negligible role: returns to M consist of a diametrical trajectory segment (that 
can be arbitrary long) and a single bounce on one of the almost flat arcs (cf. IICZ05I ). 
This allows us to compare directly the derivatives of the return maps Too and Tg in such 
points. Throughout, we use the local orthogonal section to describe tangent maps (cf. 



Appendix Bi. Consider x E M with return time n{x) — N +1 performing diametrical 



quasi-integrable motion, and let us denote the derivatives of the return maps of the c = °o 
and c = c cases by = Dv(7L) and ~ Dx{fi), respectively. Then D^r — Bt^ - A and 
Dm — Bn ■ A, where the matrix Bn is the tangent map corresponding to consecutive 
diametrical bounces on the bottom arc, while the matrices A and A are the tangent maps 
corresponding to the single bounce on the flat arc of the c = °o case and on the almost 
flat arc of the c = c case, respectively. Then we have: 

• Dn is strongly hyperbolic, expanding unstable vectors by at least a factor ki ■ N for 
some numerical constant A:i > (see IICZ05II ). 
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All elements of the matrix Bn have absolute value less than k2 ■ N for some 
numerical constant k2 > 0. (This can easily be checked by direct computation based 



on equation (B.l i.) 

• Choosing c big enough, the difference of the matrices A and A can be made 
arbitrarily small: that is, for any e > there exist some finite c' such that whenever 
c > c' we have | |A — A| | < £ (which again follows from direct computation). 
The above three facts imply that, choosing c big enough, expands unstable vectors at 
least by a factor k^, N for some positive constant ki,. In particular, at points that give rise 
to diametrical motion, uniform hyperbolicity of the return map persists for finite, large 
enough c. 

(iii) Thus we only need to show that the third type of quasi-integrable motion, sliding 
trajectories, that is, points with ^ ~ or (j? ~ tt, does not spoil the picture. 

For brevity we introduce the coordinate y/ where = 9 if 9 ~ 0, and xjf = n — (p if 
<p ~ TT. The arguments above can be summarized as follows: fix the parameters {b,c) where 
b < 1 and c > c{b). Then there exists \j/t, c such that T is uniformly hyperbolic for any x E M 
for which \j/ > \j/b,c- Furthermore, for any fixed b < 1, yf[,,c — ^ as c ^> 00. 



Analysis of sliding trajectories. To complete our heuristic argument, we show that for c 
large enough, trajectories leave the sliding region of the phase space (the region where ~ 0) 
within uniformly bounded time. 

Note that \j/ remains constant as long as the sliding trajectory collides with the same arc. 
Thus what we need to understand is how the value of y/ changes when the trajectory bypasses 
a corner point (and thus switches to a neighboring arc). There are two types of transitions to 
be considered: switching from the bottom (or the top) arc to an almost-flat arc (an /? — > 
transition), or from an almost flat arc to the bottom (or the top) arc (an Roo ~> R transition). 

Consider such a transition and denote the value of the collision angle before and after the 
corner point by y/ and y/, respectively. To describe configurations, it is more convenient to 
introduce one more quantity, j3, the angular distance of the last bounce on the first arc from 
the corner point (the distance in terms of the coordinate k divided by the radius of the circular 
arc). We have < /3 < 2y/. First let us consider an Rao ^ R transition. Using elementary 
geometry - sine theorems for the two triangles which are delimited by the radial segments 
connecting the center of the arc with the location of the last (first) bounce, the radial segment 
pointing to the corner point, and the trajectory - we find 

-cos v/'+ 1 1 — - I cos(j3 -y/) = COSY- 
c \ cj 

Now, as Y (^nd thus /3) is small, we approximate the cosines by second-order Taylor 
polynomials, and as c is large, we neglect corrections proportional to ^, to obtain 

V/'2«i/2^cj3(2v/-j3). 

This means that y' ^ cy> Si^Roo^ R transitions, with high probability - unless the last bounce 
on the almost-flat arc, or the first bounce on the bottom arc is extremely close to the corner 
point. 

On the other hand, for 7? — > transitions a completely analogous argument yields 
which means that y' ^ Y with high probability in this case. 
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As 7? — > /?oo and Roo ^ R transitions alternate, for a typical sliding trajectory the value of 
y/ exhibits a big kick at every second transition, and changes moderately at the intermediate 
occasions. This results in escaping from the sliding region after a couple of transitions. 

These phenomena have been tested by convincing simulations, which can be reproduced 
by the reader using the program available at BHall . Note that on the phase portraits produced 
by this program, k is displayed as the configurational (horizontal) coordinate instead of j3 
- this should be kept in mind when comparing the simulated trajectories with the above 
calculations. 



4.2. Lyapunov-weighted dynamics 



Further support for Conjecture |3.1| is provided by numerically searching the two-dimensional 
parameter space for pairs {b,c) where no elliptic islands are found, since the system exhibits 
a dichotomy: the system is (presumably) ergodic if there are no elliptic islands, whereas the 
presence of elliptic islands - even small ones - implies that the system is not ergodic. 

The simplest method to search for elliptic islands consists of sampling from a grid of 
initial conditions in phase space and estimating the Lyapunov exponent or an equivalent 
indicator I SESF04 I for each one. However, this is unreliable and time-consuming, since the 
number of initial conditions needed to locate an island is inversely proportional to the area of 
the island in phase space. 

A more powerful approach to locate regions of phase space with regular (non-chaotic) 
dynamics is the Lyapunov-weighted dynamics method (LWD), introduced in ref. OTK07I . The 
idea of this method is to evolve a population of walkers under a modified version of the 
dynamics, chosen so that the cloud of walkers spontaneously concentrates in regions of phase 
space which have a large (respectively small) Lyapunov exponent, according to whether a 
parameter a of the method is positive (respectively negative) fTKOTl. 

To do so, each walker follows the underlying deterministic dynamics of the system under 
study, but perturbed by a small random noise of strength e. Each walker also carries a tangent 
vector, which evolves under the tangent dynamics of the map. The local stretching factor 
of this tangent vector is calculated, as in the standard calculation of Lyapunov exponents, 
and walkers are killed or copied ("cloned") at a rate which is proportional to (a times) the 
stretching of their tangent vector IITK07L This process can be shown to lead to the desired 
effect of the cloud of walkers "highlighting" the regions of chaotic or regular behaviour, 
depending on the value of a chosen IITK07I . Details of the application of the LWD method 
to billiard models will be discussed in a forthcoming publication fHTSV 

Figure [8] illustrates the results of applying the LWD method to a generalised squash 
billiard, searching for both regular (a = — 1) and chaotic (a = +1) dynamics. WTien a = —1, 
the walkers concentrate in one of the elliptic islands of the phase space, whereas when a=l, 
the walkers instead stay in part of a chaotic component of the phase space, as shown in figure 
8(b) I In this case, there is in fact a partial barrier to transport in the phase space, visible around 



cos(^) ~ 0.7 in figure[8] The LWD dynamics with a = +1 nonetheless only sees part of this 
chaotic component. 

Our approach is thus to exhaustively explore the {b,c) parameter space, looking for 
elliptic islands by using a modified version of the LWD method IIHTSI with a < 0. As 
discussed above, in this case the walkers indeed have a tendency to concentrate in the part 



of phase space which is "most stable". Figure 9(a) shows the results of applying LWD to a 



squash billiard with parameters b = 0.1 and c = 4. For these parameter values, there is an 
elliptic periodic trajectory of period 6, around which there is an elliptic island, but there is 
also a set of parabolic periodic orbits, corresponding to diametrical bouncing, since the lower 



Chaos and stability in a two-parameter family of convex billiard tables 



14 



1,0 



0.5 



9- 

0.0 

o 
o 



-0.5 



-1.0 











0.0 0.2 0.4 0.6 0.8 1.0 
k 

(a) a = -l 




(b) a = +l 



Figure 8. Results of applying LWD to a squash billiard with parameters h = 0.6 and c = 1.6, 
searching for both (a) regular and (b) chaotic regions. The thin (blue) points represent the 
evolution of randomly chosen initial conditions that evolve with the normal dynamics of the 
squash billiard, thus giving the standard representation of the phase space of the system. The 
(red) darker points show the superposition of the phase space locations during the last 100 
collisions of 1000 walkers evolved under LWD for 10000 collisions. 
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Figure 9. Results of applying LWD with a = — 1 to two squash billiards. The superposition 
of the last 100 collisions of 10000 walkers evolved under LWD for 10000 collisions is shown 
in (red) darker points on top of the phase space. 



circular arc is longer than a semi-circle, as discussed in section |Z2| We see that when we apply 
LWD, the walkers concentrate exclusively around the elliptic island, with none trapped near 
the parabolic orbits, as shown in figure 9(a) This preference for the more stable islands is 
independent of both the size and the period of the island, as indicated in figure 9(b) where 
there are two islands of period 18, as well as parabolic bouncing orbits, for c — 6. 

If, however, we now increase c further, to c = 13, then all of the walkers concentrate 
around the parabolic orbits, as shown in figure 10(a) This strongly suggests that there are no 
longer any elliptic islands in phase space, which would be more stable, and that the system is 
thus ergodic. 

Applying LWD with a > 0, that is, searching for the chaotic region, instead shows a 
complicated fractal structure, as shown in figure 10(b) This corresponds to the region of 
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phase space where the rate of expansion is maximised. This region appears to fill phase space 
in a non-uniform way, avoiding the parabolic orbits. For comparison, figure 1 1 shows the 
results of applying LWD to a system which is known to be ergodic, namely the Bunimovich 
stadium. Similar behaviour is found as for the squash billiard, which reinforces the idea that 
the squash is ergodic. 

Nonetheless, it is, of course, still possible that small elliptic islands in phase space still 
remain in the case of the squash billiard, but whose area is below the threshold which may be 
detected using the LWD method. This threshold is the consequence of the noisy nature of the 
method, which we intend to explore in more detail in future work IHTSII . 
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Figure 10. Results of applying LWD to a squash billiard with parameters i> = 0.1 and c = 13. 
The supeiposition of the last 100 collisions of 1000 walkers evolved under LWD for 10000 
collisions is shown in (red) darker points on top of the phase space. 



Figure 12 shows a sketch of the complete two-dimensional {b,c) parameter space for 
c < 50. Here we depict the qualitative behaviour found with LWD, i.e. whether or not we 
find elliptic islands for each given pair {b,c). This picture provides strong evidence for our 



Conjecture 3. 1 confirming that no islands are found by the method for c above a certain value 
which depends on b. The disappearance and reappearance of islands observed in the figure 
as c is increased for certain values of b is presumably related to the interplay of geometric 



destruction and complex series of period-doubling bifurcations [Mac93 1 (see also section 3.3 i. 



For even larger values of c, we have also taken random values of c and b and we still 



observe the behaviour shown in figure 10 where the walkers concentrate around the parabolic 
orbits. As far as the LWD method allows us to exclude the existence of elliptic islands, we 
thus obtain further evidence pointing towards our Conjecture |3. 1| that the system is ergodic 
for c greater than some critical value c{b). 



5. Conclusions 



We have studied the dynamics of a two-parameter class of generalised squash billiards, which 
interpolates between completely integrable and completely chaotic dynamics. We have shown 
that the dynamical properties in the two-dimensional parameter plane are rather rich, involving 
a mixture of period-doubling behaviour reminiscent of smooth dynamics, and destruction of 
orbits caused by collisions with corners of the billiard table. 

We have conjectured, based on heuristic arguments which we hope can be made rigorous, 
and extensive numerical simulations with both standard and Lyapunov- weighted methods, that 
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Figure 11. Results of applying LWD to an ergodic Bunimovich stadium. The superposition 
of the last 100 collisions of 1000 walkers evolved under LWD for 10000 collisions is shown 
in (red) darker points on top of the phase space. 




parameter b 

Figure 12. Parameter space of the generalised squash billiard. LWD with 1000 walkers 
evolving during 2000 collisions is applied to find regular regions (a = —1). Pairs {b.c) that 
have elliptic islands found in this way are shown with (red) diamonds; pairs that do not have 
them are shown with (blue) points. The black continuous curve corresponds to co{b). 



all elliptic periodic points and their associated islands disappear for tables which are close 
enough to skewed stadia, thus giving a previously unknown class of ergodic convex billiards. 
It remains to characterize the parameter space in more detail, and prove the conjecture. 
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Figure Al. Geometry of generalised squash billiards. The trapezoid is shown in (blue) dashed 
lines for a generic value of h. The different circular arcs are shown in different grey scales 
(colours). 



Appendix A. Construction of generalised squasli billiards 

In this appendix, we sketch the geometric construction of generalised squash billiards. For 
fixed values of the parameters b and c, we must determine the radii and centers of the circular 
arcs that make up the table so that they satisfy the conditions of having common tangents at 



their points of intersection. To determine these quantities, we use the notation of Figure Al 
Fixing the parameter b determines the length of the top of the trapezoid and the points 

P\, P2, P3 and P4. Since the table is symmetric about the vertical axis, it is only necessary to 

satisfy the tangent conditions for three arcs, since the right and left arcs have the same radii 

and their centers are reflected in the vertical axis. 

The upper circular arc has radius r and center Cr, the lower arc has radius R and center 

Cr, and the right arc has radius Roo and center Cro^. Denote by M the midpoint of the segment 

joining Pi and P2, and by ^ a vector perpendicular to this segment. 

By symmetry, C, is on the y-axis, taking a suitable Cartesian coordinate system. Since 

II-P2 ^ C,-|| — I", we have 



Cr=[0,P2y-^r^~Pij. (A.l) 

This allows us to calculate the tangent to the upper circular arc T2 at P2, obtaining T2 = 
{Cr — P2)_L, where denotes a vector perpendicular to a given vector v. 

Let us denote, furthermore, by n2 the unit vector perpendicular to 72. Since Cr^ ~ 
P2 + sn2, with 5 e M, we have T2 ■ Cr^ — T2 P2- Moreover, since Cr^ is on the line through M 
with direction q, we find 

T2-P2-T2M 

Cr =M + toq, where t2 = — — = = . (A.2) 

T2-q 

With this we can calculate the tangent to the lower arc at the point Pi, obtaining 
ri = (C«^-Pi)±. 

Finally, Ti Cr = T\ ■ Py and Cr is on the vertical axis, so 

Cr = (0,fi), where fi = Afr-^- (A.3) 
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These equations give us Cr and Cr^ as a function of r, since T2 is a function of r and Ty 
is function of Cr^, which is also a function of r. Substituting these equations in the definition 
of the parameter c, we obtain 

c^^^^^, (A.4) 
R WCr-PiW ' 

giving an implicit equation for r in terms of b and c: 

Q^f{r,b,c)^\\CR^~P,\\~c\\CR~P2\\ . (A.5) 

This equation may be solved numerically, for example by bisection, to find the value of r 
corresponding to given values of b and c. This can then be substituted back into equations 



(A. 2 1 and (A. 3 1 to get Cr^ and Cr, and with this we obtain the radii of the other arcs; 

R^ = \\Cr^-Px\\; (A.6) 
R^\\Cr-P2\\. (A.7) 

Appendix B. Numerical methods 

For the two-parameter generalised squashes, to simulate the dynamics it is sufficient to find 
the intersection of the line representing the straight trajectory with the different circular arcs, 
which may be done by standard methods | Gas9 8|. 

A critical role in the dynamics is played by the periodic orbits, which reveal much 
dynamical information. Such orbits may be found by observing the dynamics for several 
thousand iterates, and searching for iterates whose coordinates lie in a small neighborhood of 
the initial condition. 

Once a periodic orbit is found, its stability properties can be calculated using the tangent 
map of the dynamics. To do so, we use a common convention: instead of the phase space 
coordinates, the (outgoing) local orthogonal section of the bilhard map is used IICM06I . The 
flow coordinates just after collision are taken as the location r € Q and the velocity v € S'; 
perturbations, up to linear order, are described by the coordinates {dr,dv), where both dr and 
dv are perpendicular to the velocity v. We fix a phase point x = {k,(p) EM, and denote its 
image under the dynamics as x' = Tx — [k' , <p') e M. We denote by K' the curvature of Y at 
the point k' , and by T the length of the free flight (the distance between k and k'). Then in 
the above-mentioned (outgoing) local orthogonal section coordinates, the image of a tangent 
vector {dr,dv) at x is the tangent vector {d/ ^dv') at x' , given by IICM06II 

^ d( '^L \ = \ IK' , , 'iK'r I ( ). (B.l) 




dv' J \ dv 

The tangent map for a higher iterate, and in particular for the first-return map corresponding to 
a periodic orbit, can be calculated as the product of several such matrices. Note that in the local 
orthogonal section coordinates the billiard map is area preserving, i.e. it has determinant 1. 
Thus the eigenvalue spectrum of D is characterized by the trace s of D: the orbit is hyperbolic, 
parabolic or elliptic, if \s\ > 2, \s\ — 2 or \s\ < 2, respectively. Lyapunov exponents of orbits 
can be also estimated, as the exponential growth rate of the trace of the matrix D„ giving the 
tangent map of the «th iterate. 
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